Revised Interoceptive Accuracy Scale (IAS-R)

Analysis

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)

Study 1 - Reanalysis of Murphy et al. (2020)

Reanalyze their data to confirm the factor structure of the IAS.

Data Preparation

Code
# https://osf.io/3m5nh
df1 <- haven::read_sav("../data/Murphy2020/Study 1.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

# https://osf.io/3m5nh
df2 <- haven::read_sav("../data/Murphy2020/Study 6 IAS.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other"))))

# Campos (2022) - Study 1 (https://osf.io/j6ef3)
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

# https://osf.io/3eztd
df4 <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) # No tickle because same Chinese character

# https://osf.io/3eztd
df5 <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) # No tickle because same Chinese character

df6 <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  )

df7 <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3
  ) |>
  filter(!is.na(Urinate))

Participants

  • Sample 1: Data from Murphy’s (2020) study 1, downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).
  • Sample 2: Data from Murphy’s (2020) study 6, downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).
  • Sample 3: Data from Campos’ (2022) study 1, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
  • Sample 4: Data from Lin’s (2023) study 1 and 3, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
  • Sample 5: Data from Lin’s (2023) study 2, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
  • Sample 6: New dataset from Makowski (under preparation), included 251 participants (Mean age = 37.5, SD = 13.2, range: [17, 76]; Gender: 53.0% women, 42.6% men, 2.79% non-binary, 1.59% missing; Education: Bachelor, 39.84%; Doctorate, 4.78%; High School, 31.47%; Master, 19.12%; missing, 0.40%; Other, 4.38%).
  • Sample 7: New dataset from Makowski (under preparation), included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).

Total N = 3743.

Descriptive

Distribution

Code
vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

dens1 <- estimate_density(select(df1, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1")
dens2 <- estimate_density(select(df2, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(setdiff(vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7 <- estimate_density(select(df7, all_of(setdiff(vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 7")

rbind(dens1, dens2, dens3, dens4, dens5, dens6, dens7) |>
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line() +
  theme_minimal() +
  labs(title = "Distribution of Responses", x = "Response", y = "Density", color = "Item") +
  facet_wrap(~Sample, scales = "free", nrow = 4)

Code
data1 <- normalize(select(df1, all_of(dens1$Parameter)), range = c(1, 5))
data2 <- normalize(select(df2, all_of(dens2$Parameter)), range = c(1, 5))
data3 <- normalize(select(df3, all_of(dens3$Parameter)), range = c(1, 5))
data4 <- normalize(select(df4, all_of(dens4$Parameter)), range = c(1, 5))
data5 <- normalize(select(df5, all_of(dens5$Parameter)), range = c(1, 5))
data6 <- select(df6, all_of(dens6$Parameter))
data7 <- select(df7, all_of(dens7$Parameter))

data_all <- rbind(
  data1, data2, data3,
  mutate(data4, Tickle = NA), mutate(data5, Tickle = NA),
  data6, mutate(data7, Taste = NA, Cough = NA, Blood_Sugar = NA)
)

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )
}

make_correlation(data_all)

Code
make_correlation(data1)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data4)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7)

Graph Analysis

Exploratory Graph Analysis (EGA) is a recently developed framework for psychometric assessment, that can be used to estimate the number of dimensions in questionnaire data using network estimation methods and community detection algorithms, and assess the stability of dimensions and items using bootstrapping. See https://r-ega.net/articles/ega.html for details.

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

Code
uva <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.306

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Hungry  Thirsty 0.222
 Urinate Defecate 0.219
  Sneeze    Cough 0.204
Code
uva$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva <- EGAnet::UVA(data = data1, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva <- EGAnet::UVA(data = data4, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data5, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva$keep_remove
NULL
Code
uva <- EGAnet::UVA(data = data7, cut.off = 0.3)
uva
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva$keep_remove
NULL

Networks

Code
make_egas <- function(data) {
  egas <- list()
  for (model in c("glasso", "TMFG")) {
    for (algo in c("walktrap", "louvain")) {
      for (type in c("ega", "ega.fit")) { # "hierega"
        if (type == "ega.fit" & algo == "louvain") next # Too slow
        egas[[paste0(model, "_", algo, "_", type)]] <- EGAnet::bootEGA(
          data = data,
          seed = 123,
          model = model,
          algorithm = algo,
          EGA.type = type,
          type = "resampling",
          plot.itemStability = FALSE,
          verbose = FALSE
        )
      }
    }
  }

  p <- EGAnet::compare.EGA.plots(
    egas$glasso_walktrap_ega, egas$glasso_walktrap_ega.fit,
    egas$glasso_louvain_ega, egas$TMFG_louvain_ega,
    egas$TMFG_walktrap_ega, egas$TMFG_walktrap_ega.fit,
    labels = c(
      "glasso_walktrap_ega", "glasso_walktrap_ega.fit",
      "glasso_louvain_ega", "TMFG_louvain_ega",
      "TMFG_walktrap_ega", "TMFG_walktrap_ega.fit"
    ),
    rows = 3,
    plot.all = FALSE
  )$all
  list(egas = egas, p = p)
}


egas_all <- make_egas(data_all)
egas_all$p

Code
egas1 <- make_egas(data1)
egas1$p

Code
egas2 <- make_egas(data2)
egas2$p

Code
egas3 <- make_egas(data3)
egas3$p

Code
egas4 <- make_egas(data4)
egas4$p

Code
egas5 <- make_egas(data5)
egas5$p

Code
egas6 <- make_egas(data6)
egas6$p

Code
egas7 <- make_egas(data7)
egas7$p

Structure Stability

Scores

Code
ega <- egas_all$egas$glasso_louvain_ega$EGA
plot(ega)

Code
ega_scores_1 <- EGAnet::net.scores(data1, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_2 <- EGAnet::net.scores(data2, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_3 <- EGAnet::net.scores(data3, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_4 <- EGAnet::net.scores(data4, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_5 <- EGAnet::net.scores(data5, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_6 <- EGAnet::net.scores(data6, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

ega_scores_7 <- EGAnet::net.scores(data7, ega)$scores$std.scores |>
  as.data.frame() |>
  setNames(c("EGA1", "EGA2", "EGA3", "EGA4"))

Factor Analysis

How Many Factors

Code
n <- parameters::n_factors(data_all, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4, 5), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Scree (SE) Scree_SE
5 Optimal coordinates Scree
5 Parallel analysis Scree
5 Kaiser criterion Scree
Code
n <- parameters::n_factors(data1, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data2, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
Code
n <- parameters::n_factors(data3, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
4 BIC BIC
Code
n <- parameters::n_factors(data4, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data5, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 BIC BIC
Code
n <- parameters::n_factors(data6, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Bentler Bentler
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 Scree (SE) Scree_SE
Code
n <- parameters::n_factors(data7, n_max = 12)

plot(n)

Code
knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ], row.names = FALSE)
n_Factors Method Family
1 Acceleration factor Scree
1 Scree (R2) Scree_SE
1 VSS complexity 1 VSS
1 Velicer’s MAP Velicers_MAP
4 beta Multiple_regression
4 Optimal coordinates Scree
4 Parallel analysis Scree
4 Kaiser criterion Scree
4 BIC BIC
4 BIC (adjusted) BIC

EFA

Code
efa5 <- parameters::factor_analysis(data_all, n = 5, rotation = "oblimin", sort = TRUE)
Loading required namespace: GPArotation
Code
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR3 MR2 MR1 MR5 MR4 Complexity Uniqueness
Cough 0.72 -5.89e-04 -0.06 -5.43e-03 0.03 1.02 0.50
Sneeze 0.67 -0.03 0.04 0.06 -0.05 1.04 0.53
Burp 0.55 0.13 0.06 -0.03 -0.02 1.15 0.59
Wind 0.48 0.11 -9.06e-03 -1.90e-03 0.06 1.15 0.68
Vomit 0.35 0.03 0.24 0.14 -0.09 2.32 0.68
Itch 0.05 0.71 -0.10 -0.06 0.12 1.12 0.46
Tickle 0.08 0.68 -0.03 0.06 4.30e-03 1.05 0.45
Bruise 8.13e-04 0.45 0.22 0.06 -0.08 1.56 0.68
Blood_Sugar 0.06 0.36 0.14 0.04 -0.06 1.43 0.79
Affective_touch -0.05 0.36 0.34 0.18 -0.16 2.91 0.64
Muscles 0.03 0.25 0.21 0.17 0.07 2.95 0.72
Taste 0.16 0.18 0.17 0.13 0.05 3.96 0.75
Urinate 0.10 -0.08 0.55 0.01 0.20 1.39 0.55
Defecate 0.22 -0.04 0.46 -0.01 0.15 1.69 0.60
Pain -0.01 0.22 0.28 0.11 0.23 3.22 0.65
Heart 1.28e-03 -0.03 -0.01 0.66 -0.01 1.01 0.59
Breathing 0.07 7.82e-03 -0.08 0.56 0.14 1.20 0.61
Sex_arousal 0.10 0.07 0.18 0.25 0.12 2.87 0.73
Thirsty -1.95e-04 0.05 0.11 6.26e-03 0.62 1.07 0.53
Hungry 0.01 0.04 -5.98e-03 0.12 0.60 1.08 0.55
Temp 0.11 0.12 0.14 0.15 0.27 3.12 0.68

The 5 latent factors (oblimin rotation) accounted for 38.36% of the total variance of the original data (MR3 = 10.15%, MR2 = 9.46%, MR1 = 6.49%, MR5 = 6.20%, MR4 = 6.05%).

Code
efa4 <- parameters::factor_analysis(data2, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Code
efa4 <- parameters::factor_analysis(data2, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR3 MR4 Complexity Uniqueness
Urinate 0.69 -0.18 0.13 -0.04 1.22 0.54
Defecate 0.68 -0.13 0.04 0.06 1.10 0.55
Breathing 0.56 0.03 0.06 7.96e-03 1.03 0.63
Hungry 0.55 0.15 -7.15e-03 0.02 1.16 0.60
Thirsty 0.52 0.05 0.14 -0.01 1.16 0.62
Sex_arousal 0.50 0.14 -0.04 0.20 1.49 0.56
Temp 0.46 0.23 0.16 -0.01 1.77 0.52
Pain 0.46 0.43 -0.05 -0.03 2.03 0.49
Taste 0.40 0.21 -0.07 0.12 1.77 0.69
Heart 0.40 0.06 -0.04 0.15 1.36 0.76
Muscles 0.35 0.24 0.07 0.09 2.07 0.65
Vomit 0.32 0.06 0.23 0.11 2.20 0.68
Tickle 0.02 0.69 0.01 0.04 1.01 0.47
Itch -0.02 0.67 0.18 0.03 1.16 0.41
Bruise 0.07 0.44 0.15 0.06 1.32 0.65
Blood_Sugar -0.08 0.38 0.16 0.09 1.61 0.77
Affective_touch 0.38 0.38 -0.13 0.05 2.26 0.63
Sneeze 0.05 0.05 0.71 0.01 1.02 0.41
Cough 0.04 0.02 0.71 0.04 1.01 0.43
Wind 7.52e-03 -0.04 -5.18e-03 0.94 1.00 0.15
Burp -0.01 0.27 0.17 0.47 1.91 0.48

The 4 latent factors (oblimin rotation) accounted for 44.25% of the total variance of the original data (MR1 = 17.71%, MR2 = 11.72%, MR3 = 7.65%, MR4 = 7.17%).

Code
efa4 <- parameters::factor_analysis(data3, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR2 MR4 Complexity Uniqueness
Cough 0.88 -0.08 -0.04 0.11 1.05 0.29
Sneeze 0.74 0.05 1.42e-03 -0.03 1.01 0.41
Vomit 0.55 0.12 0.04 -0.09 1.17 0.58
Burp 0.55 0.08 0.22 -0.13 1.47 0.48
Wind 0.53 0.10 0.21 -0.25 1.84 0.48
Temp 0.46 0.34 -0.06 0.09 1.97 0.52
Taste 0.29 0.18 0.10 0.12 2.37 0.74
Hungry -8.85e-03 0.70 -0.01 -0.09 1.03 0.54
Urinate 0.08 0.66 -0.12 0.11 1.15 0.52
Thirsty -0.08 0.64 -0.03 0.14 1.14 0.62
Defecate 0.10 0.63 -3.50e-05 -0.07 1.08 0.53
Breathing -0.10 0.49 0.28 -0.07 1.76 0.66
Sex_arousal 0.17 0.44 0.17 -0.22 2.20 0.58
Pain 0.21 0.41 0.08 0.18 2.00 0.58
Heart 0.05 0.28 0.24 0.08 2.23 0.77
Muscles 0.24 0.27 0.20 0.19 3.66 0.61
Tickle 0.02 -0.05 0.84 -0.04 1.01 0.31
Itch 0.03 0.02 0.74 0.11 1.05 0.39
Affective_touch 0.08 0.23 0.36 0.12 2.08 0.67
Bruise 0.17 0.07 0.29 0.45 2.08 0.55
Blood_Sugar 4.11e-03 0.11 0.35 0.43 2.07 0.60

The 4 latent factors (oblimin rotation) accounted for 45.69% of the total variance of the original data (MR1 = 15.50%, MR3 = 15.46%, MR2 = 11.02%, MR4 = 3.71%).

Code
efa5 <- parameters::factor_analysis(data4, n = 5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR5 MR4 MR2 Complexity Uniqueness
Cough 0.67 -0.07 8.90e-03 0.06 0.03 1.04 0.54
Sneeze 0.63 -0.04 0.11 0.05 -0.08 1.11 0.54
Burp 0.58 0.18 -0.03 -0.09 0.03 1.25 0.59
Wind 0.54 0.08 -0.06 -0.03 0.11 1.17 0.67
Vomit 0.33 0.06 0.21 0.15 -0.19 3.02 0.70
Taste 0.24 0.19 -4.60e-03 0.06 0.18 3.00 0.77
Bruise -0.03 0.63 0.05 0.05 -0.03 1.03 0.58
Affective_touch 0.07 0.53 0.03 -0.01 -0.04 1.05 0.67
Blood_Sugar 0.07 0.43 0.15 0.10 -0.18 1.82 0.67
Itch 0.11 0.41 -0.01 0.05 0.15 1.45 0.71
Muscles 0.07 0.40 -0.05 0.02 0.19 1.56 0.75
Urinate -0.04 0.06 0.64 6.75e-03 0.06 1.04 0.55
Defecate 0.11 0.02 0.59 -0.03 0.02 1.08 0.57
Heart -0.04 0.04 -0.03 0.66 -0.04 1.02 0.59
Breathing 0.09 -2.40e-03 -0.02 0.51 0.12 1.17 0.65
Sex_arousal 0.07 0.18 0.07 0.24 0.13 2.86 0.78
Thirsty 0.05 -0.05 0.24 0.03 0.49 1.51 0.60
Hungry 0.09 -0.10 0.15 0.21 0.38 2.25 0.67
Temp 0.09 0.17 7.51e-03 0.11 0.36 1.81 0.72
Pain -3.56e-03 0.30 0.14 0.04 0.33 2.38 0.67

The 5 latent factors (oblimin rotation) accounted for 35.08% of the total variance of the original data (MR1 = 10.08%, MR3 = 8.61%, MR5 = 5.94%, MR4 = 5.33%, MR2 = 5.12%).

Code
efa5 <- parameters::factor_analysis(data5, n = 5, rotation = "oblimin", sort = TRUE)
plot(efa5)

Code
display(efa5)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR4 MR2 MR3 MR5 Complexity Uniqueness
Burp 0.69 0.06 0.03 -0.10 0.03 1.06 0.50
Wind 0.62 0.01 0.03 -0.02 -0.08 1.04 0.64
Cough 0.54 -0.06 0.01 0.17 0.11 1.31 0.58
Sneeze 0.45 0.02 -0.07 0.21 0.17 1.78 0.61
Bruise -0.04 0.60 -0.01 0.02 0.06 1.03 0.63
Affective_touch 0.08 0.53 -0.04 -0.08 0.06 1.12 0.67
Blood_Sugar 0.01 0.48 -0.09 0.08 0.19 1.45 0.68
Itch 0.18 0.38 0.08 0.12 -0.02 1.77 0.69
Muscles 0.09 0.38 0.19 3.00e-03 -0.07 1.68 0.77
Thirsty 0.03 -0.17 0.57 -0.01 0.22 1.50 0.55
Temp 2.56e-03 0.18 0.49 0.03 -0.08 1.34 0.71
Hungry 0.02 -0.02 0.48 0.16 0.09 1.32 0.65
Taste 0.19 0.22 0.32 0.05 -0.12 2.84 0.73
Pain 0.10 0.24 0.28 2.96e-03 0.03 2.26 0.77
Breathing 0.04 -0.03 0.03 0.67 -0.07 1.03 0.54
Heart -0.08 0.05 5.68e-03 0.58 0.06 1.08 0.66
Sex_arousal 0.11 0.14 0.23 0.26 -0.06 3.08 0.76
Urinate -0.04 0.11 0.11 -0.01 0.61 1.14 0.54
Defecate 0.14 0.03 0.03 -0.02 0.60 1.11 0.53
Vomit 0.15 0.15 -0.04 0.22 0.23 3.56 0.74

The 5 latent factors (oblimin rotation) accounted for 35.25% of the total variance of the original data (MR1 = 8.99%, MR4 = 8.07%, MR2 = 6.47%, MR3 = 5.97%, MR5 = 5.75%).

Code
efa4 <- parameters::factor_analysis(data6, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR2 MR4 MR3 Complexity Uniqueness
Hungry 0.69 0.22 -0.10 -0.15 1.36 0.47
Urinate 0.68 -0.13 -0.05 0.13 1.16 0.55
Thirsty 0.57 0.18 0.13 -0.08 1.34 0.47
Temp 0.56 -0.08 0.28 0.04 1.50 0.48
Pain 0.53 0.03 0.25 0.10 1.50 0.44
Defecate 0.49 -0.05 -0.05 0.31 1.75 0.60
Sex_arousal 0.46 0.07 -0.01 0.24 1.58 0.61
Affective_touch 0.43 0.18 0.18 -0.06 1.77 0.61
Muscles 0.36 0.02 0.26 0.26 2.71 0.51
Tickle 0.07 0.66 0.04 -2.12e-03 1.03 0.49
Itch 5.18e-04 0.65 -0.05 0.12 1.09 0.56
Bruise 0.02 0.59 0.11 -0.11 1.14 0.59
Burp -0.02 0.44 0.04 0.32 1.86 0.61
Taste 0.13 0.33 0.17 0.32 2.82 0.50
Blood_Sugar 0.14 0.25 0.12 -0.09 2.44 0.85
Breathing -0.03 0.03 0.81 2.25e-03 1.01 0.35
Heart 0.12 0.02 0.50 -8.27e-04 1.12 0.66
Vomit 0.19 0.14 0.27 0.16 3.06 0.66
Sneeze 0.03 0.28 0.12 0.45 1.84 0.54
Wind 0.18 0.14 0.15 0.42 1.88 0.54
Cough 0.14 0.21 0.22 0.25 3.53 0.62

The 4 latent factors (oblimin rotation) accounted for 44.31% of the total variance of the original data (MR1 = 16.65%, MR2 = 11.34%, MR4 = 9.45%, MR3 = 6.86%).

Code
efa4 <- parameters::factor_analysis(data7, n = 4, rotation = "oblimin", sort = TRUE)
plot(efa4)

Code
display(efa4)
Rotated loadings from Factor Analysis (oblimin-rotation)
Variable MR1 MR3 MR4 MR2 Complexity Uniqueness
Hungry 0.62 -0.04 0.15 0.04 1.13 0.51
Sex_arousal 0.58 0.06 0.01 0.03 1.03 0.60
Pain 0.55 0.04 0.09 -8.09e-03 1.06 0.61
Thirsty 0.55 0.09 0.07 0.03 1.10 0.58
Temp 0.36 0.11 0.19 0.09 1.85 0.65
Affective_touch 0.28 -0.05 0.13 0.26 2.51 0.75
Defecate 0.06 0.67 -7.30e-03 0.02 1.02 0.49
Sneeze -0.18 0.53 0.29 0.06 1.85 0.59
Vomit 0.13 0.53 -0.05 0.09 1.20 0.61
Urinate 0.40 0.47 -0.12 -0.02 2.09 0.52
Burp 2.70e-03 0.46 0.17 0.08 1.34 0.64
Breathing 0.02 0.02 0.71 0.01 1.00 0.45
Heart 0.17 7.44e-03 0.59 0.01 1.16 0.51
Muscles 0.15 0.07 0.47 0.07 1.31 0.60
Itch -5.06e-03 -0.02 -0.02 0.78 1.00 0.43
Tickle -0.01 0.02 -0.01 0.73 1.00 0.47
Bruise 8.63e-03 0.04 0.06 0.38 1.08 0.81
Wind 0.10 0.08 0.10 0.26 1.82 0.82

The 4 latent factors (oblimin rotation) accounted for 40.83% of the total variance of the original data (MR1 = 12.34%, MR3 = 10.09%, MR4 = 9.29%, MR2 = 9.12%).

CFA

Structure

Code
# Ambiguous: Temp + Vomit + Affective_touch + Sex_arousal + Taste
# Discard: Tickle + Blood_Sugar

m1 <- "
Interoception =~ Hungry + Thirsty + Urinate + Defecate + Itch + Bruise + Muscles + Pain + Breathing + Heart + Cough + Sneeze + Wind + Burp
"

m4 <- "
Sustenance =~ Hungry + Thirsty + Urinate + Defecate + Muscles + Pain
Skin =~ Itch + Bruise
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m5 <- "
Sustenance =~ Hungry + Thirsty + Pain
UrinateDefecate =~ Urinate + Defecate
Skin =~ Itch + Bruise + Muscles
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"


m6 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
Expulsion =~ Cough + Sneeze + Wind + Burp
"

m7 <- "
HungryThirsty =~ Hungry + Thirsty
UrinateDefecate =~ Urinate + Defecate
ItchBruise =~ Itch + Bruise
MusclesPain =~ Muscles + Pain
HeartBreathing =~ Breathing + Heart
CoughSneeze =~ Cough + Sneeze
WindBurp =~ Wind + Burp
"

m1_all <- lavaan::cfa(m1, data = data_all)
m4_all <- lavaan::cfa(m4, data = data_all)
m5_all <- lavaan::cfa(m5, data = data_all)
m6_all <- lavaan::cfa(m6, data = data_all)
m7_all <- lavaan::cfa(m7, data = data_all)

anova(m1_all, m4_all, m5_all, m6_all, m7_all) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
m7_all -14218.06 -13919.70 239.58 56
m6_all -14036.69 -13774.87 432.95 62 < .001
m5_all -13745.07 -13513.69 734.56 67 < .001
m4_all -13389.39 -13182.37 1098.24 71 < .001
m1_all -12372.54 -12202.05 2127.10 77 < .001

Anova Table (Type 1 tests)

Code
parameters::parameters(m7_all, standardize=TRUE) |>
  display()
# Loading
Link Coefficient SE 95% CI z p
HungryThirsty =~ Hungry 0.66 0.02 (0.63, 0.70) 40.45 < .001
HungryThirsty =~ Thirsty 0.69 0.02 (0.66, 0.73) 42.07 < .001
UrinateDefecate =~ Urinate 0.69 0.02 (0.66, 0.72) 45.07 < .001
UrinateDefecate =~ Defecate 0.71 0.02 (0.68, 0.74) 46.26 < .001
ItchBruise =~ Itch 0.59 0.02 (0.55, 0.62) 32.53 < .001
ItchBruise =~ Bruise 0.56 0.02 (0.52, 0.59) 30.99 < .001
MusclesPain =~ Muscles 0.58 0.02 (0.55, 0.61) 36.37 < .001
MusclesPain =~ Pain 0.68 0.02 (0.64, 0.71) 42.39 < .001
HeartBreathing =~ Breathing 0.70 0.02 (0.66, 0.74) 34.75 < .001
HeartBreathing =~ Heart 0.56 0.02 (0.52, 0.60) 29.78 < .001
CoughSneeze =~ Cough 0.73 0.01 (0.70, 0.75) 55.74 < .001
CoughSneeze =~ Sneeze 0.71 0.01 (0.69, 0.74) 54.21 < .001
WindBurp =~ Wind 0.67 0.01 (0.64, 0.69) 46.98 < .001
WindBurp =~ Burp 0.74 0.01 (0.71, 0.77) 53.17 < .001
# Correlation
Link Coefficient SE 95% CI z p
HungryThirsty ~~ UrinateDefecate 0.59 0.02 (0.55, 0.64) 26.16 < .001
HungryThirsty ~~ ItchBruise 0.49 0.03 (0.43, 0.54) 16.82 < .001
HungryThirsty ~~ MusclesPain 0.61 0.02 (0.56, 0.66) 24.64 < .001
HungryThirsty ~~ HeartBreathing 0.51 0.03 (0.46, 0.57) 19.82 < .001
HungryThirsty ~~ CoughSneeze 0.47 0.02 (0.43, 0.52) 20.06 < .001
HungryThirsty ~~ WindBurp 0.44 0.02 (0.39, 0.49) 17.95 < .001
UrinateDefecate ~~ ItchBruise 0.45 0.03 (0.39, 0.51) 15.70 < .001
UrinateDefecate ~~ MusclesPain 0.63 0.02 (0.58, 0.67) 26.14 < .001
UrinateDefecate ~~ HeartBreathing 0.48 0.03 (0.43, 0.53) 18.57 < .001
UrinateDefecate ~~ CoughSneeze 0.55 0.02 (0.51, 0.59) 25.07 < .001
UrinateDefecate ~~ WindBurp 0.50 0.02 (0.45, 0.54) 21.56 < .001
ItchBruise ~~ MusclesPain 0.82 0.03 (0.77, 0.88) 29.25 < .001
ItchBruise ~~ HeartBreathing 0.52 0.03 (0.46, 0.58) 17.22 < .001
ItchBruise ~~ CoughSneeze 0.64 0.03 (0.59, 0.69) 24.48 < .001
ItchBruise ~~ WindBurp 0.67 0.03 (0.62, 0.73) 25.64 < .001
MusclesPain ~~ HeartBreathing 0.56 0.03 (0.50, 0.61) 20.38 < .001
MusclesPain ~~ CoughSneeze 0.57 0.02 (0.52, 0.62) 23.69 < .001
MusclesPain ~~ WindBurp 0.59 0.02 (0.54, 0.63) 24.12 < .001
HeartBreathing ~~ CoughSneeze 0.51 0.02 (0.46, 0.55) 20.42 < .001
HeartBreathing ~~ WindBurp 0.42 0.03 (0.37, 0.47) 16.21 < .001
CoughSneeze ~~ WindBurp 0.76 0.02 (0.73, 0.80) 41.66 < .001
Code
anova(
  lavaan::cfa(m1, data = data1),
  lavaan::cfa(m4, data = data1),
  lavaan::cfa(m5, data = data1),
  lavaan::cfa(m6, data = data1),
  lavaan::cfa(m7, data = data1)
  ) |>
  parameters::parameters() |>
  display()
Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
                is not positive definite;
                use lavInspect(fit, "cov.lv") to investigate.
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data1) -1950.88 -1749.41 88.22 56
lavaan::cfa(m6, data = data1) -1924.80 -1748.01 126.29 62 < .001
lavaan::cfa(m5, data = data1) -1859.71 -1703.48 201.38 67 < .001
lavaan::cfa(m4, data = data1) -1834.28 -1694.49 234.82 71 < .001
lavaan::cfa(m1, data = data1) -1660.86 -1545.74 420.24 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m1, data = data2),
  lavaan::cfa(m4, data = data2),
  lavaan::cfa(m5, data = data2),
  lavaan::cfa(m6, data = data2),
  lavaan::cfa(m7, data = data2)
  ) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data2) -1898.38 -1705.96 66.50 56
lavaan::cfa(m6, data = data2) -1822.70 -1653.84 154.19 62 < .001
lavaan::cfa(m5, data = data2) -1800.91 -1651.68 185.98 67 < .001
lavaan::cfa(m4, data = data2) -1754.66 -1621.14 240.22 71 < .001
lavaan::cfa(m1, data = data2) -1632.46 -1522.51 374.42 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m1, data = data3),
  lavaan::cfa(m4, data = data3),
  lavaan::cfa(m5, data = data3),
  lavaan::cfa(m6, data = data3),
  lavaan::cfa(m7, data = data3)
  ) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data3) -3587.97 -3380.01 114.03 56
lavaan::cfa(m6, data = data3) -3507.56 -3325.06 206.45 62 < .001
lavaan::cfa(m5, data = data3) -3452.59 -3291.31 271.41 67 < .001
lavaan::cfa(m4, data = data3) -3418.95 -3274.65 313.05 71 < .001
lavaan::cfa(m1, data = data3) -3143.19 -3024.35 600.82 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m1, data = data4),
  lavaan::cfa(m4, data = data4),
  lavaan::cfa(m5, data = data4),
  lavaan::cfa(m6, data = data4),
  lavaan::cfa(m7, data = data4)
  ) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data4) -5136.94 -4888.93 115.89 56
lavaan::cfa(m6, data = data4) -5103.03 -4885.39 161.81 62 < .001
lavaan::cfa(m5, data = data4) -5041.07 -4848.74 233.76 67 < .001
lavaan::cfa(m4, data = data4) -4930.21 -4758.13 352.62 71 < .001
lavaan::cfa(m1, data = data4) -4591.03 -4449.31 703.81 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m1, data = data5),
  lavaan::cfa(m4, data = data5),
  lavaan::cfa(m5, data = data5),
  lavaan::cfa(m6, data = data5),
  lavaan::cfa(m7, data = data5)
  ) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data5) -1737.85 -1531.33 83.61 56
lavaan::cfa(m6, data = data5) -1711.33 -1530.10 122.14 62 < .001
lavaan::cfa(m5, data = data5) -1683.25 -1523.10 160.21 67 < .001
lavaan::cfa(m4, data = data5) -1629.89 -1486.59 221.58 71 < .001
lavaan::cfa(m1, data = data5) -1485.55 -1367.54 377.91 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m1, data = data6),
  lavaan::cfa(m4, data = data6),
  lavaan::cfa(m5, data = data6),
  lavaan::cfa(m6, data = data6),
  lavaan::cfa(m7, data = data6)
  ) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data6) -1437.92 -1265.17 91.47 56
lavaan::cfa(m6, data = data6) -1444.97 -1293.38 96.42 62 0.551
lavaan::cfa(m5, data = data6) -1412.50 -1278.53 138.89 67 < .001
lavaan::cfa(m4, data = data6) -1419.48 -1299.62 139.91 71 0.907
lavaan::cfa(m1, data = data6) -1342.29 -1243.58 229.10 77 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(str_remove(m5, fixed("Cough + ")), data = data7),
  lavaan::cfa(str_remove(m1, fixed("Cough + ")), data = data7),
  lavaan::cfa(str_remove(m4, fixed("Cough + ")), data = data7),
  lavaan::cfa(str_remove(m6, fixed("Cough + ")), data = data7),
  lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data7)
  ) |>
  parameters::parameters() |>
  display()
Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
                is not positive definite;
                use lavInspect(fit, "cov.lv") to investigate.
Warning in lavTestLRT(object = object, ..., model.names = NAMES): lavaan
WARNING: some models are based on a different set of observed variables
Parameter AIC BIC Chi2 df p
c(“lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")),”, ” data = data7)“) -1606.97 -1443.78 82.59 39
lavaan::cfa(str_remove(m6, fixed(“Cough +”)), data = data7) -1825.70 -1654.15 109.52 50 0.005
lavaan::cfa(str_remove(m5, fixed(“Cough +”)), data = data7) -1821.75 -1671.12 123.48 55 0.016
lavaan::cfa(str_remove(m4, fixed(“Cough +”)), data = data7) -1752.45 -1618.56 200.78 59 < .001
lavaan::cfa(str_remove(m1, fixed(“Cough +”)), data = data7) -1665.13 -1556.35 300.09 65 < .001

Anova Table (Type 1 tests)

Higher Order Factors

Code
m7h1 <- paste0(m7, "
              Interoception =~ HungryThirsty + UrinateDefecate + ItchBruise + MusclesPain + HeartBreathing + CoughSneeze + WindBurp
              ")

m7h3 <- paste0(m7, "
              Valenced =~ MusclesPain + HeartBreathing + ItchBruise
              Expulsion =~ CoughSneeze + WindBurp
              Homeostasis =~ HungryThirsty + UrinateDefecate
              ")

m7h1_all <- lavaan::cfa(m7h1, data = data_all)
m7h3_all <- lavaan::cfa(m7h3, data = data_all)
anova(m7_all, m7h1_all, m7h3_all) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
m7_all -14218.06 -13919.70 239.58 56
m7h3_all -14104.40 -13873.02 375.24 67 < .001
m7h1_all -13903.22 -13690.11 582.42 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data1), 
  lavaan::cfa(m7h1, data = data1), 
  lavaan::cfa(m7h3, data = data1)) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data1) -1950.88 -1749.41 88.22 56
lavaan::cfa(m7h3, data = data1) -1903.59 -1747.35 157.51 67 < .001
lavaan::cfa(m7h1, data = data1) -1873.38 -1729.48 193.71 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data2), 
  lavaan::cfa(m7h1, data = data2), 
  lavaan::cfa(m7h3, data = data2)) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data2) -1898.38 -1705.96 66.50 56
lavaan::cfa(m7h3, data = data2) -1855.12 -1705.90 131.76 67 < .001
lavaan::cfa(m7h1, data = data2) -1841.29 -1703.85 151.59 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data3), 
  lavaan::cfa(m7h1, data = data3), 
  lavaan::cfa(m7h3, data = data3)) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data3) -3587.97 -3380.01 114.03 56
lavaan::cfa(m7h3, data = data3) -3580.71 -3419.43 143.30 67 0.002
lavaan::cfa(m7h1, data = data3) -3496.70 -3348.16 233.30 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data4), 
  lavaan::cfa(m7h1, data = data4), 
  lavaan::cfa(m7h3, data = data4)) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data4) -5136.94 -4888.93 115.89 56
lavaan::cfa(m7h3, data = data4) -5103.39 -4911.06 171.44 67 < .001
lavaan::cfa(m7h1, data = data4) -5007.06 -4829.92 273.77 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data5), 
  lavaan::cfa(m7h1, data = data5), 
  lavaan::cfa(m7h3, data = data5)) |>
  parameters::parameters() |>
  display()
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data5) -1737.85 -1531.33 83.61 56
lavaan::cfa(m7h3, data = data5) -1727.18 -1567.03 116.28 67 < .001
lavaan::cfa(m7h1, data = data5) -1709.86 -1562.35 139.61 70 < .001

Anova Table (Type 1 tests)

Code
anova(
  lavaan::cfa(m7, data = data6), 
  lavaan::cfa(m7h1, data = data6), 
  lavaan::cfa(m7h3, data = data6)) |>
  parameters::parameters() |>
  display()
Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
                is not positive definite;
                use lavInspect(fit, "cov.lv") to investigate.
Parameter AIC BIC Chi2 df p
lavaan::cfa(m7, data = data6) -1437.92 -1265.17 91.47 56
lavaan::cfa(m7h3, data = data6) -1438.54 -1304.58 112.85 67 0.030
lavaan::cfa(m7h1, data = data6) -1430.27 -1306.88 127.12 70 0.003

Anova Table (Type 1 tests)

None of the model can be properly estimated (variances are negative).

Code
NA
[1] NA
Code
# anova(
#   lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data7), 
#   lavaan::cfa(str_remove(str_remove(m7h1, fixed("\nCoughSneeze =~ Cough + Sneeze")), fixed("CoughSneeze +")), data = data7), 
#   lavaan::cfa(str_remove(str_remove(m7h3, fixed("\nCoughSneeze =~ Cough + Sneeze")), fixed("CoughSneeze +")), data = data7)) |>
#   parameters::parameters() |>
#   display()

Model Performance

Code
rbind(
  performance::performance(m7_all, metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "All"),
  performance::performance(update(m7_all, data = data1), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 1"),
  performance::performance(update(m7_all, data = data2), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 2"),
  performance::performance(update(m7_all, data = data3), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 3"),
  performance::performance(update(m7_all, data = data4), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 4"),
  performance::performance(update(m7_all, data = data5), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 5"),
  performance::performance(update(m7_all, data = data6), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 6"),
  performance::performance(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data7), metrics = c("Chi2", "RMSEA", "CFI", "SRMR")) |>
    mutate(Sample = "Sample 7")
) |>
  display()
Chi2 RMSEA CFI SRMR Sample
239.58 0.03 0.98 0.02 All
88.22 0.04 0.98 0.03 Sample 1
66.50 0.02 0.99 0.03 Sample 2
114.03 0.04 0.98 0.03 Sample 3
115.89 0.03 0.98 0.02 Sample 4
83.61 0.03 0.98 0.03 Sample 5
91.47 0.05 0.97 0.04 Sample 6
82.59 0.05 0.97 0.03 Sample 7

No evidence for higher order factors.

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

Code
scores <- list(
  sample1 = cbind(ega_scores_1, as.data.frame(predict(lavaan::cfa(m7, data = data1)))) |>
    mutate(Sample = "Sample 1"),
  sample2 = cbind(ega_scores_2, as.data.frame(predict(lavaan::cfa(m7, data = data2)))) |>
    mutate(Sample = "Sample 2"),
  sample3 = cbind(ega_scores_3, as.data.frame(predict(lavaan::cfa(m7, data = data3)))) |>
    mutate(Sample = "Sample 3"),
  sample4 = cbind(ega_scores_4, as.data.frame(predict(lavaan::cfa(m7, data = data4)))) |>
    mutate(Sample = "Sample 4"),
  sample5 = cbind(ega_scores_5, as.data.frame(predict(lavaan::cfa(m7, data = data5)))) |>
    mutate(Sample = "Sample 5"),
  sample6 = cbind(ega_scores_6, as.data.frame(predict(lavaan::cfa(m7, data = data6)))) |>
    mutate(Sample = "Sample 6"),
  sample7 = cbind(ega_scores_7, as.data.frame(predict(lavaan::cfa(str_remove(m7, fixed("\nCoughSneeze =~ Cough + Sneeze")), data = data7)))) |>
    mutate(Sample = "Sample 7", CoughSneeze = NA)
)
  
save(scores, file = "../data/scores.RData")

Study 2 - Correlates

Code
load("../data/scores.RData")

df1 <- cbind(df1, scores$sample1)
df2 <- cbind(df2, scores$sample2)
df3 <- cbind(df3, scores$sample3)
df4 <- cbind(df4, scores$sample4)
df5 <- cbind(df5, scores$sample5)
df6 <- cbind(df6, scores$sample6)
df7 <- cbind(df7, scores$sample7)

vars_intero <- names(select(scores$sample1, -Sample))
vars_cfa <- names(select(scores$sample1, -Sample, -starts_with("EGA")))
vars_cfa7 <- vars_cfa[vars_cfa != "CoughSneeze"]

# names(df1)
# names(df2)
# names(df3)  # TAS, BPQ
# names(df4)  # BPQ
# names(df5)  # PHQ, TAS
# names(df6)  # BDI, STAI
# names(df7)  # MAIA, PHQ4, PI, LIE, IPIP, PID, SPQ

EGA and CFA

Code
df_scores <- rbind(
  select(df1, all_of(vars_intero)),
  select(df2, all_of(vars_intero)),
  select(df3, all_of(vars_intero)),
  select(df4, all_of(vars_intero)),
  select(df5, all_of(vars_intero)),
  select(df6, all_of(vars_intero)),
  select(df7, all_of(vars_intero))
)

correlation::correlation(df_scores, redundant = TRUE) |>
    correlation::cor_sort() |>
    correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
    mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Parameter1, y = Parameter2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = val), size = 3) +
    labs(title = "Correlation Matrix") +
    scale_fill_metro_c(limit = c(0, 1), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )

Demographics

Code
make_diff <- function(df, vars, component=NA) {
  dat <- data.frame()
  for (resp in vars) {
    f <- as.formula(paste0(resp, "~ Gender"))
    d <- filter(df, Gender %in% c("Male", "Female")) |> 
      mutate(Gender = fct_relevel(Gender, "Male", "Female"))
    dat <- effectsize::cohens_d(f, data=d) |> 
      mutate(var_Correlate = "Gender", var_Interoception = resp, 
             p=t.test(f, data=d)$p.val) |> 
      as.data.frame() |> 
      rbind(dat)
  }
    
  dat$r <- effectsize::d_to_r(dat$Cohens_d, n1=sum(df$Gender=="Male"), n2=sum(df$Gender=="Female"))
  dat$CI_low <- effectsize::d_to_r(dat$CI_low, n1=sum(df$Gender=="Male"), n2=sum(df$Gender=="Female"))
  dat$CI_high <- effectsize::d_to_r(dat$CI_high, n1=sum(df$Gender=="Male"), n2=sum(df$Gender=="Female"))
  dat$Sample <- unique(df$Sample)
  dat$Component <- component
  select(dat, -Cohens_d)
}

data <- rbind(
  make_diff(df1, vars_cfa, component="Demographic"),
  make_diff(df2, vars_cfa, component="Demographic"),
  make_diff(df3, vars_cfa, component="Demographic"),
  make_diff(df4, vars_cfa, component="Demographic"),
  make_diff(df5, vars_cfa, component="Demographic"),
  make_diff(df6, vars_cfa, component="Demographic"),
  make_diff(df7, vars_cfa7, component="Demographic")
) |> 
  mutate(var_Correlate = "Gender (Male vs. Female)")
Code
make_cor <- function(df, vars, with = "Age", component=NA) {
  dat <- data.frame() 
  for (v in vars) {
    dat <- correlation::cor_test(df, v, with) |> 
      rename(var_Interoception = Parameter1, var_Correlate = Parameter2) |> 
      as.data.frame() |> 
      rbind(dat)
  }
    
  dat$Sample <- unique(df$Sample)
  dat$Component <- component
  select(dat, -df_error, -t, -Method, -n_Obs)
}

data <- rbind(
  data,
  make_cor(df1, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df2, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df3, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df4, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df5, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df6, vars_cfa, with = "Age", component="Demographic"),
  make_cor(df7, vars_cfa7, with = "Age", component="Demographic")
) 

Body

Code
# Reverse (TODO: Ana confirm why)
df3$BPQ_BodyAwareness <- -1*normalize(df3$BPQ_BodyAwareness)
df3$BPQ_AutonomicReactivity <- -1*normalize(df3$BPQ_AutonomicReactivity)

data <- rbind(
  data,
  make_cor(df3, vars_cfa, with = "BPQ_BodyAwareness", component="Body") |> 
    mutate(var_Correlate = "Body Awareness"),
  make_cor(df3, vars_cfa, with = "BPQ_AutonomicReactivity", component="Body") |> 
    mutate(var_Correlate = "Autonomic Reactivity"),
  make_cor(df4, vars_cfa, with = "BPQ_A", component="Body") |> 
    mutate(var_Correlate = "Body Awareness")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "MAIA_AttentionRegulation", component="MAIA")  |> 
    mutate(var_Correlate = "Attention Regulation"),
  make_cor(df7, vars_cfa7, with = "MAIA_BodyListening", component="MAIA") |> 
    mutate(var_Correlate = "Body Listening"),
  make_cor(df7, vars_cfa7, with = "MAIA_EmotionalAwareness", component="MAIA") |> 
    mutate(var_Correlate = "Emotional Awareness"),
  make_cor(df7, vars_cfa7, with = "MAIA_NotDistracting", component="MAIA") |> 
    mutate(var_Correlate = "Not Distracting"),
  make_cor(df7, vars_cfa7, with = "MAIA_Noticing", component="MAIA") |> 
    mutate(var_Correlate = "Noticing"),
  make_cor(df7, vars_cfa7, with = "MAIA_NotWorrying", component="MAIA") |> 
    mutate(var_Correlate = "Not Worrying"),
  make_cor(df7, vars_cfa7, with = "MAIA_SelfRegulation", component="MAIA") |> 
    mutate(var_Correlate = "Self Regulation"),
  make_cor(df7, vars_cfa7, with = "MAIA_Trusting", component="MAIA") |> 
    mutate(var_Correlate = "Trusting")
)
Code
# names(df3)[str_detect(names(df3), "TAS")]
# names(df5)[str_starts(names(df5), "TAS")]

data <- rbind(
  data,
  make_cor(df3, vars_cfa, with = "TAS_IdentifyingFeelings", component="Body") |> 
    mutate(var_Correlate = "Difficulties Identifying Feelings"),
  make_cor(df3, vars_cfa, with = "TAS_DescribingFeelings", component="Body") |> 
    mutate(var_Correlate = "Difficulties Describing Feelings"),
  make_cor(df3, vars_cfa, with = "TAS_ExternallyOrientedThinking", component="Body") |> 
    mutate(var_Correlate = "External Thinking"),
  make_cor(df5, vars_cfa, with = "TAS_DIF", component="Body") |> 
    mutate(var_Correlate = "Difficulties Identifying Feelings"),
  make_cor(df5, vars_cfa, with = "TAS_DDF", component="Body") |> 
    mutate(var_Correlate = "Difficulties Describing Feelings")
)

Mood

Code
data <- rbind(
  data,
  make_cor(df5, vars_cfa, with = "PHQ15_sum", component="Mood") |> 
    mutate(var_Correlate = "Somatic Concerns (PHQ-15)"),
  make_cor(df5, vars_cfa, with = "PHQ9_sum", component="Mood")|> 
    mutate(var_Correlate = "Depression (PHQ-9)"),
  make_cor(df7, vars_cfa7, with = "PHQ4_Depression", component="Mood") |> 
    mutate(var_Correlate = "Depression (PHQ-4)"),
  make_cor(df7, vars_cfa7, with = "PHQ4_Anxiety", component="Mood") |> 
    mutate(var_Correlate = "Anxiety (PHQ-4)")
)
Code
data <- rbind(
  data,
  make_cor(df6, vars_cfa, with = "BDI2_Total", component="Mood") |> 
    mutate(var_Correlate = "Depression (BDI-II)")
)
Code
data <- rbind(
  data,
  make_cor(df6, vars_cfa, with = "STAI5_General", component="Mood") |> 
    mutate(var_Correlate = "Anxiety (STAI-5)")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "BPD", component="Mood") |> 
    mutate(var_Correlate = "Boderline Personality")
)

Psychopathology

Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "PID5_Antagonism", component="Maladaptive") |> 
    mutate(var_Correlate = "Antagonism"),
  make_cor(df7, vars_cfa7, with = "PID5_Detachment", component="Maladaptive") |> 
    mutate(var_Correlate = "Detachment"),
  make_cor(df7, vars_cfa7, with = "PID5_Disinhibition", component="Maladaptive") |> 
    mutate(var_Correlate = "Disinhibition"),
  make_cor(df7, vars_cfa7, with = "PID5_NegativeAffect", component="Maladaptive") |> 
    mutate(var_Correlate = "Negative Affect"),
  make_cor(df7, vars_cfa7, with = "PID5_Psychoticism", component="Maladaptive") |> 
    mutate(var_Correlate = "Psychoticism")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "SPQ_ConstrictedAffect", component="Schizotypic") |> 
    mutate(var_Correlate = "Constricted Affect"),
  make_cor(df7, vars_cfa7, with = "SPQ_Eccentric", component="Schizotypic") |> 
    mutate(var_Correlate = "Eccentric"),
  make_cor(df7, vars_cfa7, with = "SPQ_MagicalThinking", component="Schizotypic") |> 
    mutate(var_Correlate = "Magical Thinking"),
  make_cor(df7, vars_cfa7, with = "SPQ_NoCloseFriends", component="Schizotypic") |> 
    mutate(var_Correlate = "No Close Friends"),
  make_cor(df7, vars_cfa7, with = "SPQ_OddSpeech", component="Schizotypic") |> 
    mutate(var_Correlate = "Odd Speech"),
  make_cor(df7, vars_cfa7, with = "SPQ_Reference", component="Schizotypic") |> 
    mutate(var_Correlate = "Reference"),
  make_cor(df7, vars_cfa7, with = "SPQ_SocialAnxiety", component="Schizotypic") |> 
    mutate(var_Correlate = "Social Anxiety"),
  make_cor(df7, vars_cfa7, with = "SPQ_Suspiciousness", component="Schizotypic") |> 
    mutate(var_Correlate = "Suspiciousness"),
  make_cor(df7, vars_cfa7, with = "SPQ_UnusualPerceptions", component="Schizotypic") |> 
    mutate(var_Correlate = "Unusual Perceptions")
)
Code
#  names(select(df7, starts_with("ASQ")))
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "ASQ_Imagination", component="Autistic") |> 
    mutate(var_Correlate = "Imagination"),
  make_cor(df7, vars_cfa7, with = "ASQ_LackSocialSkills", component="Autistic") |>
    mutate(var_Correlate = "Lack of Social Skills"),
  make_cor(df7, vars_cfa7, with = "ASQ_LowAttentionalSwitching", component="Autistic") |>
    mutate(var_Correlate = "Low Attentional Switching"),
  make_cor(df7, vars_cfa7, with = "ASQ_Patterns", component="Autistic") |>
    mutate(var_Correlate = "Patterns and Numbers"),
  make_cor(df7, vars_cfa7, with = "ASQ_Routine", component="Autistic") |>
    mutate(var_Correlate = "Routines")
)

Personality

Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "IPIP6_Agreeableness", component="Personality") |> 
    mutate(var_Correlate = "Agreeableness"),
  make_cor(df7, vars_cfa7, with = "IPIP6_Conscientiousness", component="Personality") |>
    mutate(var_Correlate = "Conscientiousness"),
  make_cor(df7, vars_cfa7, with = "IPIP6_Neuroticism", component="Personality") |>
    mutate(var_Correlate = "Neuroticism"),
  make_cor(df7, vars_cfa7, with = "IPIP6_Openness", component="Personality") |>
    mutate(var_Correlate = "Openness"),
  make_cor(df7, vars_cfa7, with = "IPIP6_Extraversion", component="Personality") |>
    mutate(var_Correlate = "Extraversion"),
  make_cor(df7, vars_cfa7, with = "IPIP6_HonestyHumility", component="Personality") |>
    mutate(var_Correlate = "Honesty-Humility")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "PI_Alive", component="World Beliefs") |> 
    mutate(var_Correlate = "Alive"),
  make_cor(df7, vars_cfa7, with = "PI_Changing", component="World Beliefs") |>
    mutate(var_Correlate = "Changing"),
  make_cor(df7, vars_cfa7, with = "PI_Enticing", component="World Beliefs") |>
    mutate(var_Correlate = "Enticing"),
  make_cor(df7, vars_cfa7, with = "PI_Good", component="World Beliefs") |>
    mutate(var_Correlate = "Good"),
  make_cor(df7, vars_cfa7, with = "PI_Hierarchical", component="World Beliefs") |>
    mutate(var_Correlate = "Hierarchical"),
  make_cor(df7, vars_cfa7, with = "PI_Safe", component="World Beliefs") |>
    mutate(var_Correlate = "Safe")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "GCBS_Extraterrestrial", component="Conspiracy Beliefs") |> 
    mutate(var_Correlate = "Extraterrestrial"),
  make_cor(df7, vars_cfa7, with = "GCBS_GlobalConspiracies", component="Conspiracy Beliefs") |>
    mutate(var_Correlate = "Global Conspiracies"),
  make_cor(df7, vars_cfa7, with = "GCBS_GovernmentMalfeasance", component="Conspiracy Beliefs") |>
    mutate(var_Correlate = "Government Malfeasance"),
  make_cor(df7, vars_cfa7, with = "GCBS_InformationControl", component="Conspiracy Beliefs") |>
    mutate(var_Correlate = "Information Control"),
  make_cor(df7, vars_cfa7, with = "GCBS_PersonalWellbeing", component="Conspiracy Beliefs") |>
    mutate(var_Correlate = "Personal Wellbeing")
)
Code
data <- rbind(
  data,
  make_cor(df7, vars_cfa7, with = "LIE_Ability", component="Lying Profile") |> 
    mutate(var_Correlate = "Ability"),
  make_cor(df7, vars_cfa7, with = "LIE_Frequency", component="Lying Profile") |>
    mutate(var_Correlate = "Frequency"),
  make_cor(df7, vars_cfa7, with = "LIE_Contextuality", component="Lying Profile") |>
    mutate(var_Correlate = "Contextuality"),
  make_cor(df7, vars_cfa7, with = "LIE_Negativity", component="Lying Profile") |>
    mutate(var_Correlate = "Negativity")
)

Summary

Code
p <- data |>
  mutate(sig = ifelse(p < .01, "p < .01", ifelse(p < .001, "p < .001", "N.S.")),
         dir = sign(r),
         var_Interoception = str_replace(var_Interoception, "CoughSneeze", "Cough/Sneeze"),
         var_Interoception = str_replace(var_Interoception, "HeartBreathing", "Heart/Breathing"),
         var_Interoception = str_replace(var_Interoception, "HungryThirsty", "Hungry/Thirsty"),
         var_Interoception = str_replace(var_Interoception, "ItchBruise", "Itch/Bruise"),
         var_Interoception = str_replace(var_Interoception, "MusclesPain", "Muscles/Pain"),
         var_Interoception = str_replace(var_Interoception, "UrinateDefecate", "Urinate/Defecate"),
         var_Interoception = str_replace(var_Interoception, "WindBurp", "Wind/Burp"),
         var_Interoception = str_replace(var_Interoception, "Expulsion_Sudden", "Expulsion - Sudden"),
         var_Interoception = fct_relevel(var_Interoception, "Heart/Breathing", "Muscles/Pain", "Itch/Bruise", "Hungry/Thirsty", "Urinate/Defecate", "Wind/Burp", "Cough/Sneeze"),
         Component = fct_relevel(
           Component, "Demographic", "MAIA", "Body", "Mood", "Maladaptive",
           "Schizotypic", "Autistic", "Personality", "Conspiracy Beliefs", "World Beliefs", "Lying Profile")) |>
  ggplot(aes(x = var_Correlate)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  # geom_bar(aes(fill=Sample, alpha=sig, y=r, group=Sample), stat="identity", position = position_dodge2(width = 0.8)) +
  geom_linerange(aes(color=dir, alpha=sig, ymin=0, ymax=r, group=Sample), position = position_dodge2(width = 0.7), linewidth=2.5) +
  geom_linerange(aes(ymin = CI_low, ymax = CI_high, color = dir, alpha=sig, group=Sample), position = position_dodge2(width = 0.7)) +
  scale_fill_gradient(low = "#2196F3", high = "#FF5722", guide = "none") +
  scale_color_gradient(low = "#2196F3", high = "#FF5722", guide = "none") +
  scale_alpha_discrete(range = c(0.2, 0.7), guide = "none") +
  scale_y_continuous(expand = c(0, 0), breaks = c(-0.3, 0, 0.3)) +
  facet_grid(Component ~ var_Interoception, scales = "free_y", switch = "y") +
  coord_flip(ylim = c(-0.45, 0.45)) +
  theme_modern(axis.title.space = 5) +
  theme(panel.grid.major.y = element_line(),
        axis.title.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.y = element_text(size = rel(0.5)),
        axis.text.x = element_text(size = rel(0.8)),
        plot.title = element_text(face = "bold", hjust = 0),
        strip.placement = "outside",
        strip.background.y = element_rect(fill = "#E0E0E0", color = "white"),
        strip.background.x = element_rect(fill = "#F8BBD0", color = "white")) 
Warning: Using alpha for a discrete variable is not advised.
Code
ggsave("figures/Figure2.png", p, width=21*0.6, height=29.7*0.45, dpi=200, bg="white")
p

Discussion

Benefits of the IAS: - Straightforward and sensation-centered items

Recommendations: - Remove Itch (redundant + issue in Chinese) - Use analog scales

Limitations: - Not much clear theorethical or empirical structure (small grouping of items) - Limited variability (clear mode at 4/5) - Ambiguous items which grouping depends on the context (and its awareness) - E.g., heart beating fast, vomit when nauseaous - Few items for some modalities (e.g., 1 for heart) - Positive phrasing: benefits but also might exacerbate positivity bias (and thus unidimensionality)

Need for context-specific items (cross-modal when possible, i.E., cardioception, respiroception, etc.).

New Scale: Multimodal Interoceptive Sensitivity Scale (MISS)